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1 Introduction 

EEG/MEG noise has an unequal variance and is 
correlated, both in space and in time. Noise variance 
may differ greatly between samples or sensors, and 
correlations between samples or sensors can be very 
high [1-4]. 

If these noise characteristics are neglected, then an 
EEG/MEG source analysis will yield unreliable 
results [e.g. 5, 6], First, source parameter estimates 
will be inefficient. That is, their standard errors will 
be too high. Second, the estimated covariance matrix 
of the parameter estimates will be inaccurate. In 
general it will give a too optimistic impression of 
precision. Third, goodness of fit measures will be 
unreliable, which may result in over- or 
undermodeling of the data. 

For these reasons, it is very beneficial to incorporate 
the spatiotemporal noise covariance in the analysis. 
Although the spatial covariance is incorporated quite 
often [6-11], the temporal covariance is disregarded 
up to now. Therefore, we are developing a method 
to incorporate the spatiotemporal noise covariance 
matrix. The essential feature of this method is that 
the estimation problem is split into two parts. First a 
model is fitted to the spatiotemporal noise 
covariance matrix. Then the source parameters are 
estimated given this noise model. 

2 Methods 

The regression model is: 

Y = FM + E (1) 

Y is the N by T matrix with EEG and/or MEG 
measurements on N sensors and T samples. F is the 
N by D matrix with unit activity of D sources. M is 
the D by T matrix with source amplitudes. E is the N 
by T matrix with noise. Vec(E), all the elements of 
E stored in a vector, is distributed with zero mean 
and NT by NT covariance matrix V. In general V is 
estimated from trial variation around the mean. 
Generalized Least Squares (GLS) estimates can be 
obtained by minimizing: 

[vec(FM) - vec(Y)] T V _1 [vec(FM) - vec(Y)] 


( 2 ) 

where T denotes transposition, and "* matrix 
inversion. Such an analysis has three disadvantages. 
First, V is singular if the number of trials does not 
exceed NT. Second, the estimate of V is imprecise, 
due to the poor ratio between the number of trials 
and the number of elements in V [11]. Third, V is 
large; therefore source estimation is computationally 
intensive. 

These problems can be circumvented if V is 
modeled by a matrix with three characteristics: a) it 
is nonsingular; b) it only depends on a few 
parameters, and thus is accurate even if only a few 
trials are measured, c) its structure simplifies source 
estimation considerably. 

2.1 Noise covariance model 

V is modeled as: 

V=U®W (3) 

That is, as the Kronecker product [e.g. 1, 12 ] of a 
temporal Tby T covariance matrix U and a spatial N 
by N covariance matrix W. In this manner, the 
number of free parameters in V is greatly reduced. A 
further reduction can be obtained by invoking 
parametric models for the temporal and spatial 
matrices. 

The temporal covariance matrix is modeled by an 
autoregressive moving average (ARMA) process of 
limited order. It has been shown that such a model 
can give an adequate description of the temporal 
noise covariance. Moreover, such a model only 
requires a few parameters, namely the sum of: a) the 
number of AR parameters, b) the number of MA 
parameters, c) one parameter for the variance of the 
moving average regressors. 

The spatial covariance matrix is allowed to be 
heteroscedastic, and the correlations are modeled as 
functions of the distance between sensors. Let a nn 
be the standard deviation of the noise on sensor n, 
and let d , be the distance between sensors n and 
n ’ then these function are [cf. 13]: 



( 8 ) 


w nn =<y nn <7 n , n ,exp(-d nn ,/a) 

w nn , =o nn o n ^xp(-d nn J a) P 

w nn < =G nn G n , rC {\ - od nn ) 

w nn , = a nn a n v ( 1 - ad nn , ) P (4) 

where a and [5 determine the size of the correlations. 
These models are chosen because noise is spatially 
heteroscedastic [3], and because correlations are 
mainly functions of the distance between sensors [1, 
3], Moreover, the number of parameters is quite low, 
namely the sum of: a) standard deviation parameters 
for every sensor, and b) one or two parameters 
modeling the correlations. 

Since EEG is measured with respect to a reference, 
the spatial covariance matrix also depends on the 
reference [3], Therefore, the aforementioned spatial 
covariance functions have to be transformed to the 
referenced case, as outlined in [3], 

2.2 Estimation of noise parameters 

Maximum likelihood estimates of the noise 
parameters are obtained by minimizing [14]: 

ln|V| - ln|V| + tr[VV“ l ]-TN (5) 

which is equivalent to minimizing [12]: 

Mn|U| + 7Tn| W| + tr[V U" '®W 1 ] ( 6 ) 

Note that there is an indeterminacy in eq. 3: 
multiplying the temporal covariance matrix by a 
constant, and dividing the spatial covariance matrix 
by the same constant yields the same results. 
Therefore, eq. 6 is minimized subject to the 
constraint that the sum of the spatial noise variances 
equals one. 

It is necessary to determine the appropriate order of 
the ARMA process, and the appropriate functions 
for the spatial correlations. The Akaike Information 
Criterion (AIC) [14] is used for this purpose. An 
efficient procedure is to determine the temporal 
model first for a few randomly chosen sensors, and 
then determine the spatial model for a few randomly 
chosen samples. 

2.3 Estimation of source parameters 

Maximum likelihood estimates of the source 
parameters can then be derived by minimizing [cf. 
15]: 


tr^U 1 [FM - Y] 'W ’[FM - Y] 


where tr denotes the trace. Note that eq. 8 can be 
much easier evaluated than eq. 2. 

3 Results 

We investigated whether the noise covariance model 
gives an adequate description of real data. The data 
pertain to noise on 64 EEG sensors referenced to the 
left mastoid, 40 samples, and 213 trials. 

3.1 Description of the data 

Noise variance. The temporal model requires that 
the noise variance is homoscedastic in time, whereas 
the spatial model allows the variances to be 
heteroscedastic in space. The variances are indeed 
temporally homogeneous, and spatially 
inhomogeneous (Fig. la). 

Spatial covariance. The Kronecker product requires 
that the spatial covariances are temporally 
homogeneous. This is indeed the case. The spatial 
model requires that the spatial correlations are 
merely a function of distance, which is indeed the 
case (Fig. lb). 

Temporal covariance. The Kronecker product 
requires that the lagged autocorrelations are spatially 
homogeneous, which is indeed the case. The 
temporal model requires that the lagged 
autocorrelations are stationary, which is also the 
case. The oscillatory pattern in the autocorrelations 
suggests that they indeed may be described by an 
ARMA model (Fig. lc). 

Spatiotemporal covariance. The Kronecker product 
implies that the lagged crosscovariance between 
sensors n and iT can be derived by multiplying the 
autocovariance function by the spatial covariance. 
The data suggest that this may indeed be appropriate 
(Fig. Id). 

3.2 Fitting the model 

As a first step, the models U and V are selected by 
analyzing three randomly chosen sensors and three 
randomly chosen samples separately 
For each sensor, the AIC indicates that an 
AR(2)MA(2) model should be selected. The 
parameter estimates are very homogeneous over 
sensors. For each sample, the exponential two 
parameter model is selected. The parameter 
estimates are very homogeneous over samples. 
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Figure 1: (a) the percentage of sensors with a 

certain average noise variance. The average noise 
variance is computed over samples, (b) average 
noise correlations as a function of distance between 
sensors. The average is computed over samples, (c) 
average autocorrelation averaged over sensors and 
lags, (d) average lagged cross correlation for 
sensors which correlate high, medium or low at zero 
lag. The average is computed over lags. 


This tentative analysis suggests that the proposed 
model may adequately describe the data. If the 
model indeed proves to be appropriate, then this will 
be very beneficial. In that case the modeled matrix is 
nonsingular, only depends on a few parameters, and 
offers the possibility to obtain source estimates in a 
computationally efficient manner. 

Future developments will include: a) comparison 
between maximum likelihood and OLS source 
estimates; b) efficient implementation using the fact 
that source strengths appear linearly; c) goodness of 
fit measure to determine the required number of 
sources; d) computation of the source parameter 
covariance matrix. 
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